##################################################################
# File:      BJPS_Replication.R
# Purporse:  Replication of main text figures and tables for 
#                 British Journal of Political Science
# Author:    Harry Oppenheimer
# Date:      02/27/25
# Last Edit:  
# Data In:   ~/Replication
# Data Out:  ~/Replication/MainText
# Prev file: None
# Status:    In Progress.
# Computer: 2.6 GHz 6-Core Intel Core i7 laptop with 32 GB of DDR4 RAM, MacOS version 13.5, R 4.2.1
##################################################################

rm(list = ls())

#Load in the libraries for the analysis
libraries <- c(  "data.table","tidyverse", "fixest", "ggpubr")
new.packages <- libraries[!(libraries %in% installed.packages()[,"Package"])]
if(length(new.packages)>0){install.packages(new.packages)} 
lapply(libraries, require, character.only = TRUE)

#Set the seed for the analysis
set.seed(02138)

#List the dictionary to make the tables later
table_dictionary = c("allied"= "Treaty",
                     "binary_pta" = "Binary PTA",
                     "wto_joint" = "Joint WTO", 
                     "time_seq" = "Month",
                     "pair" = "Dyad",
                     "paird" = "Directed Dyad",
                     "from" = "CountryA",
                     "to" = "CountryB",
                     "n_agreements_all" = "Interconnection Agreements",
                     "diff_n_agreements_all" = "Differenced Interconnection Agreements",
                     "fromt"= "CountryA * Month",
                     "tot" = "CountryB * Month",
                     "log_agreements_all" = "log(Agreements + 1)")

#Set the working directory to the replication files location
setwd("~/Dropbox/Research/Projects/Replication")
dir.create("MainText") 

data <- fread("bjps_data_all.csv")
data_soe <- fread("bjps_data_SOE.csv")


###############
#Summary statistics discussed in the subsection "Measuring Digital Interdependence"
###############

#Number of countries and dyads
length(unique(data$from))
length(unique(data$paird))

#Standard deviation, mean, largest value
sd(data$n_agreements_all)
mean(data$n_agreements_all)
data$pair[which.max(data$n_agreements_all)]

#Proportion of all 0 dyads
q <- data %>% group_by(paird) %>%
  summarize(mall = mean(n_agreements_all) > 0)
prop.table(table(q$mall))

#Proportion of 0s in dyads with at least one non 0

q2 <- data %>% filter(paird %in% q$paird[which(q$mall == T)])
prop.table(table(q2$n_agreements_all > 0))


###############
#Code to create Figure 1 in the main text 
###############

n_interconnect <- data %>%
  group_by(month) %>% summarise(n_all = sum(n_agreements_all, na.rm = T),
                                n_pair = n_distinct(paird[n_agreements_all > 0])) %>%
  mutate(pct_ch_agreements = (n_all - lag(n_all)) / lag(n_all) * 100)

mean(n_interconnect$pct_ch_agreements, na.rm = T)

p1 <- ggplot(n_interconnect, aes(x=month)) + geom_line(aes(y = n_pair), size = 1) + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank()) + ggtitle("Dyads With at Least One Agreement") + 
  scale_y_continuous(breaks=seq(2600, 10000, 400))

p2 <- ggplot(n_interconnect, aes(x=month)) + geom_line(aes(y = n_all), size = 1) + 
  theme(axis.title.x=element_blank(),
        axis.title.y=element_blank()) + ggtitle("Total Number of International Agreements") + 
  scale_y_continuous(breaks=seq(50000, 250000, 25000))
ggpubr::ggarrange(p1, p2) 
ggsave("MainText/Figure1.png", 
       width = 8, height  = 4,
       device = "png", dpi = 600)


###############
#Code for the analysis in Table 1 in the main text
###############

m1 = feglm(n_agreements_all ~ allied | paird  +  fromt  + tot,
                 data,
                 family=quasipoisson,
                 cluster = ~paird+time_seq)

m2 = feglm(n_agreements_all ~ allied  + wto_joint | paird  +  fromt  + tot,
                  data,
                     family=quasipoisson,
                     cluster = ~paird+time_seq)

m3 = feglm(n_agreements_all ~ allied  + binary_pta | paird  +  fromt  + tot,
                       data,
                       family=quasipoisson,
                       cluster = ~paird+time_seq)

m4 = feglm(n_agreements_all ~ allied  + binary_pta + wto_joint | paird  +  fromt  + tot,
                     data,
                     family=quasipoisson,
                     cluster = ~paird+time_seq)

etable(m1, m2,  m3, m4, tex = T, digits = 4, dict = table_dictionary)

###############
#Code for the analysis in Table 2 in the main text 
###############

data_subset <- data[which(data$ixp_joint_2010 == 1 | data$select2_2010 == 1),]

m1 = feglm(n_agreements_all ~ allied | paird  +  fromt  + tot,
           data_subset,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m2 = feglm(n_agreements_all ~ allied  + wto_joint | paird  +  fromt  + tot,
           data_subset,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m3 = feglm(n_agreements_all ~ allied  + binary_pta | paird  +  fromt  + tot,
           data_subset,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m4 = feglm(n_agreements_all ~ allied  + binary_pta + wto_joint | paird  +  fromt  + tot,
           data_subset,
           family=quasipoisson,
           cluster = ~paird+time_seq)

etable(m1, m2,  m3, m4, tex = T, digits = 4, dict = table_dictionary)

###############
#Code for the analysis in Table 3 in the main text 
###############

m1 = feglm(n_agreements_all ~ allied | paird  +  fromt  + tot,
           data_soe,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m2 = feglm(n_agreements_all ~ allied  + wto_joint | paird  +  fromt  + tot,
           data_soe,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m3 = feglm(n_agreements_all ~ allied  + binary_pta | paird  +  fromt  + tot,
           data_soe,
           family=quasipoisson,
           cluster = ~paird+time_seq)

m4 = feglm(n_agreements_all ~ allied  + binary_pta + wto_joint | paird  +  fromt  + tot,
           data_soe,
           family=quasipoisson,
           cluster = ~paird+time_seq)

etable(m1, m2,  m3, m4, tex = T, digits = 4, dict = table_dictionary)

###############
#Z tests in the main text 
###############

#Z test comparing coefficient in Table 1 Model 1 to Table 2 Model 1
(0.4606 - 0.4532) / (.0839^2 + .0846^2)^.5

#Z test comparing coefficient in Table 1 Model 1 to Table 3 Model 1
(0.4606 - .3950) / (.0839^2 + .1031^2)^.5



